################################################################################
## Paper:     Working the Crowd
## Authors:   P. Mongrain, N. Fréchet, B. Thompson Collart, and Y. Dufresne
## Date:      February 2025
################################################################################

################################################################################
## IMPORTANT NOTE
################################################################################

# Please run the "merge.do" *before* running the code

################################################################################
## SET WORKING DIRECTORY
################################################################################

getwd()

setwd("C:/...") #add appropriate file path

################################################################################
## LOAD PACKAGES 
################################################################################

#remotes::install_github("NightingaleHealth/ggforestplot")

library(cowplot)
library(devtools)
library(dfoptim)
library(dotwhisker)
library(dplyr) 
library(effects)
library(emmeans)
library(expss)
library(ggalt)
library(ggeasy)
library(ggeffects)
library(ggforestplot)
library(ggplot2)
library(ggpubr)
library(ggridges)
library(grid)
library(gridExtra)
library(haven)
library(interplot)
library(lattice)
library(lme4)
library(lubridate)
library(magrittr)
library(marginaleffects)
library(margins)
library(merTools)
library(prettyR)
library(psych)
library(rms)
library(scales)
library(sjmisc)
library(sjPlot)
library(srvyr)
library(tibble)
library(tidyverse)
library(visreg)

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, lme4, marginaleffects, optimx, broom.mixed)

################################################################################
## LOAD DATA
################################################################################

## IMPORT DATA FROM STATA

merge <- read_dta("merge.dta")

## CREATE SEPARATE DATAFRAME FOR EACH ELECTION

ca2011.data <- merge %>% filter(election == "ca2011")
ca2015.data <- merge %>% filter(election == "ca2015")
ca2019.data <- merge %>% filter(election == "ca2019")
on2011.data <- merge %>% filter(election == "on2011")
on2014.data <- merge %>% filter(election == "on2014")
qc2022.data <- merge %>% filter(election == "qc2022")

## DEFINE VARIABLES AS FACTOR OR NUMERIC

qc2022.data$correct_district <- as.factor(qc2022.data$correct_district)
qc2022.data$vote_district <- as.factor(qc2022.data$vote_district)
qc2022.data$pidstatus_district <- as.factor(qc2022.data$pidstatus_district)
qc2022.data$univ <- as.factor(qc2022.data$univ)
qc2022.data$male <- as.factor(qc2022.data$male)
qc2022.data$age55 <- as.factor(qc2022.data$age55)
qc2022.data$highinc <- as.factor(qc2022.data$highinc)
qc2022.data$reelected <- as.factor(qc2022.data$reelected)
qc2022.data$margin <- as.numeric(qc2022.data$margin)
qc2022.data$interest_n <- as.numeric(qc2022.data$interest_n)
qc2022.data$vote_district <- factor(qc2022.data$vote_district, levels = c(0,1),
                                    labels = c("No", "Yes"))

## CREATE MODE FUNCTION

getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

################################################################################
## RANDOM INTERCEPT
################################################################################

## Estimate model

qc2022.data = qc2022.data %>%
    mutate(margin = (margin - mean(margin)) / sd(margin)) #scaling margin

qc2022.data = qc2022.data %>%
  mutate(time = (time - mean(time)) / sd(time)) #scaling response date

pacman::p_load(lmerTest) #show statistical significance

qc2022.ri.fit <- glmer(correct_district ~ 1 + vote_district*univ + male + 
                     age55 + highinc + reelected + margin + time +
                     (1 | district_code), weight=survey_weight, 
                     data=qc2022.data, 
                     family=binomial(link="logit"))

summary(qc2022.ri.fit, correlation = FALSE)

qc2022.ri.int1 <- emmeans(qc2022.ri.fit, ~ univ | vote_district) #first difference
qc2022.ri.int2 <- pairs(pairs(emmeans::regrid(qc2022.ri.int1)), by = NULL) #second difference

qc2022.ri.int.em1 <- emmeans(qc2022.ri.fit, pairwise ~ univ | vote_district)
qc2022.ri.int.em2 <- pairs(qc2022.ri.int.em1) %>%  confint()

## Extracting random effects

ranef(qc2022.ri.fit, 
      drop = TRUE, #to have them as a vector
      condVar = FALSE)

## Plotting random effects

broom.mixed::tidy(qc2022.ri.fit, effects = "ran_vals") %>%
    ggplot(aes(x = reorder(level, -estimate), y = estimate)) +
    geom_hline(yintercept = 0, color = "grey70") +
    geom_point(color = "cyan4") +
    geom_errorbar(aes(ymin = estimate - 1.96*std.error, ymax = estimate + 1.96*std.error), 
                  width = 0, color = "cyan4") +
    labs(x = "District", y = "RE estimate") +
    scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
    theme_classic() +
    coord_flip()

################################################################################
## RANDOM SLOPE
################################################################################

qc2022.rs.fit.bob <- glmer(correct_district ~ 1 + vote_district*univ + male + 
                    age55 + highinc + margin + time + reelected + 
                    (1 + vote_district*univ | district_code), weight=survey_weight,
                    control = glmerControl(optimizer = "bobyqa"),
                    data=qc2022.data, 
                    family=binomial(link="logit"))

summary(qc2022.rs.fit.bob, correlation = FALSE)

qc2022.rs.int1 <- emmeans(qc2022.rs.fit.bob, ~ univ | vote_district) #first difference
qc2022.rs.int2 <- pairs(pairs(emmeans::regrid(qc2022.rs.int1)), by = NULL) #second difference

qc2022.rs.int.em1 <- emmeans(qc2022.rs.fit.bob, pairwise ~ univ | vote_district)
qc2022.rs.int.em2 <- pairs(qc2022.rs.int.em1) %>% confint()

################################################################################
## PREDICTED PROBABILITIES
################################################################################

qc2022.data <- filter(qc2022.data, !is.na(correct_district))
qc2022.data <- filter(qc2022.data, !is.na(univ))
qc2022.data <- filter(qc2022.data, !is.na(vote_district))
qc2022.data <- filter(qc2022.data, !is.na(male))
qc2022.data <- filter(qc2022.data, !is.na(age55))
qc2022.data <- filter(qc2022.data, !is.na(highinc))
qc2022.data <- filter(qc2022.data, !is.na(reelected))
qc2022.data <- filter(qc2022.data, !is.na(margin))

qc2022.rs.preds.1 <- with(qc2022.data, data.frame(district_code=unique(district_code), 
                                                  univ="0", 
                                                  interest_n=0, 
                                                  vote_district="No",
                                                  male=getmode(male),
                                                  age55=getmode(age55),
                                                  highinc=getmode(highinc),
                                                  reelected=getmode(reelected), 
                                                  margin=0,
                                                  time=0))

qc2022.rs.preds.2 <- with(qc2022.data, data.frame(district_code=unique(district_code), 
                                                  univ="1", 
                                                  interest_n=1, 
                                                  vote_district="No",
                                                  male=getmode(male),
                                                  age55=getmode(age55),
                                                  highinc=getmode(highinc),
                                                  reelected=getmode(reelected), 
                                                  margin=0,
                                                  time=0))

qc2022.rs.preds.3 <- with(qc2022.data, data.frame(district_code=unique(district_code), 
                                                  univ="0", 
                                                  interest_n=0, 
                                                  vote_district="Yes",
                                                  male=getmode(male),
                                                  age55=getmode(age55),
                                                  highinc=getmode(highinc),
                                                  reelected=getmode(reelected), 
                                                  margin=0,
                                                  time=0))

qc2022.rs.preds.4 <- with(qc2022.data, data.frame(district_code=unique(district_code), 
                                                  univ="1", 
                                                  interest_n=1, 
                                                  vote_district="Yes",
                                                  male=getmode(male),
                                                  age55=getmode(age55),
                                                  highinc=getmode(highinc),
                                                  reelected=getmode(reelected), 
                                                  margin=0,
                                                  time=0))

qc2022.rs.preds.1$preds = predict(qc2022.rs.fit.bob, newdata = qc2022.rs.preds.1, type = "response")
qc2022.rs.preds.2$preds = predict(qc2022.rs.fit.bob, newdata = qc2022.rs.preds.2, type = "response")
qc2022.rs.preds.3$preds = predict(qc2022.rs.fit.bob, newdata = qc2022.rs.preds.3, type = "response")
qc2022.rs.preds.4$preds = predict(qc2022.rs.fit.bob, newdata = qc2022.rs.preds.4, type = "response")

qc2022.rs.preds.1 <- qc2022.rs.preds.1 %>% mutate(pr.type = 1)
qc2022.rs.preds.1 <- qc2022.rs.preds.1 %>% arrange(preds)
qc2022.rs.preds.1$id <- seq(1, nrow(qc2022.rs.preds.1))
order.district.loser <- subset(qc2022.rs.preds.1, select = c(district_code, id))

qc2022.rs.preds.2.2 <- qc2022.rs.preds.2 %>% mutate(pr.type = 88)
qc2022.rs.preds.2 <- qc2022.rs.preds.2 %>% mutate(pr.type = 2)
qc2022.rs.preds.2 <- qc2022.rs.preds.2 %>% arrange(preds)
qc2022.rs.preds.2$id <- seq(1, nrow(qc2022.rs.preds.2))

qc2022.rs.preds.2.2 <- merge(qc2022.rs.preds.2.2, order.district.loser, by="district_code")
qc2022.rs.preds.2.2$id <- seq(1, nrow(qc2022.rs.preds.2.2))

qc2022.rs.preds.3 <- qc2022.rs.preds.3 %>% mutate(pr.type = 3)
qc2022.rs.preds.3 <- qc2022.rs.preds.3 %>% arrange(preds)
qc2022.rs.preds.3$id <- seq(1, nrow(qc2022.rs.preds.3))
order.district.ind <- subset(qc2022.rs.preds.3, select = c(district_code, id))

qc2022.rs.preds.4.2 <- qc2022.rs.preds.4 %>% mutate(pr.type = 99)
qc2022.rs.preds.4 <- qc2022.rs.preds.4 %>% mutate(pr.type = 4)
qc2022.rs.preds.4 <- qc2022.rs.preds.4 %>% arrange(preds)
qc2022.rs.preds.4$id <- seq(1, nrow(qc2022.rs.preds.4))

qc2022.rs.preds.4.2 <- merge(qc2022.rs.preds.4.2, order.district.ind, by="district_code")
qc2022.rs.preds.4.2$id <- seq(1, nrow(qc2022.rs.preds.4.2))

qc2022.rs.preds.loser <- rbind(qc2022.rs.preds.1, qc2022.rs.preds.2, qc2022.rs.preds.2.2)
qc2022.rs.preds.winner <- rbind(qc2022.rs.preds.3, qc2022.rs.preds.4, qc2022.rs.preds.4.2)

qc2022.rs.preds.loser$pr.type <- as.factor(qc2022.rs.preds.loser$pr.type)
qc2022.rs.preds.winner$pr.type <- as.factor(qc2022.rs.preds.winner$pr.type)

################################################################################
## FIGURE 3: QUEBEC 2022
################################################################################

qc2022.rs.preds.loser$title <- "(a) Losers"

qc2022.rs.preds.loser.g <- ggplot(qc2022.rs.preds.loser, aes(x = reorder(id, -id), 
                    y = preds)) + geom_point(aes(color = pr.type, alpha = pr.type), show.legend = FALSE) + 
                    labs(x = "District", y = "Predicted Probability of Correct Forecast") +
                    scale_x_discrete(guide = guide_axis(n.dodge = 2), expand = expansion(mult = c(.02, .02))) +
                    scale_y_continuous("Predicted Probability", limits=c(0,1), breaks = seq(0,1,.25)) +
                    geom_hline(yintercept = 0.5, color = "grey70", linetype = 2) +
                    scale_colour_manual(name = "University", values = c("#0072B2","#88CCEE","#88CCEE"),
                    labels = c("No", "Yes", "")) +
                    scale_alpha_manual(values = c(1, 1, .15)) + 
                    theme_bw() + coord_flip() + facet_grid(. ~ title)

qc2022.rs.preds.loser.g <- qc2022.rs.preds.loser.g + theme(panel.grid.major = element_blank(),
                                             panel.grid.minor = element_blank(),
                                             panel.background = element_blank(),
                                             axis.text.y = element_blank(), 
                                             axis.ticks.y = element_blank(),
                                             axis.title.y = element_blank(),
                                             axis.title.x = element_blank(),
                                             axis.text.x = element_text(size = 14, colour="white"),
                                             axis.ticks = element_line(color = "#00000000"),                                             
                                             strip.text = element_text(size = 15),
                                             plot.title = element_text(size = 12, hjust = 0.5))

qc2022.rs.preds.loser.g

qc2022.rs.preds.winner$title <- "(b) Winners"

qc2022.rs.preds.winner.g <- ggplot(qc2022.rs.preds.winner, aes(x = reorder(id, -id), 
                     y = preds)) + geom_point(aes(color = pr.type, alpha = pr.type), show.legend = FALSE) + 
                     labs(x = "District", y = "Predicted Probability of Correct Forecast") +
                     scale_x_discrete(guide = guide_axis(n.dodge = 2), expand = expansion(mult = c(.02, .02))) +
                     scale_y_continuous("Predicted Probability", limits=c(0,1), breaks = seq(0,1,.25)) +
                     geom_hline(yintercept = 0.5, color = "grey70", linetype = 2) +
                     scale_colour_manual(name = "University", values = c("#0072B2","#88CCEE","#88CCEE"),
                     labels = c("No", "Yes", "")) +
                     scale_alpha_manual(values = c(1, 1, .15)) + 
                     theme_bw() + coord_flip() + facet_grid(. ~ title)

qc2022.rs.preds.winner.g <- qc2022.rs.preds.winner.g + theme(panel.grid.major = element_blank(),
                                               panel.grid.minor = element_blank(),
                                               panel.background = element_blank(),
                                               axis.text.y = element_blank(), 
                                               axis.ticks.y = element_blank(),
                                               axis.title.y = element_blank(),
                                               axis.title.x = element_blank(),
                                               axis.text.x = element_text(size = 14, colour="white"),
                                               axis.ticks = element_line(color = "#00000000"),
                                               strip.text = element_text(size = 15),
                                               plot.title = element_text(size = 12, hjust = 0.5))

qc2022.rs.preds.winner.g

tiff("qc2022_preds_g.tiff", units="in", width=8, height=4, res=300)

qc2022.preds.g <- ggarrange(qc2022.rs.preds.loser.g, 
                           qc2022.rs.preds.winner.g,
                           ncol = 2, nrow = 1)

qc2022.preds.g <- annotate_figure(qc2022.preds.g,
                                 left = textGrob("Electoral District", rot = 90, 
                                                 vjust = .6, gp = gpar(col = "white", cex = 1)),
                                 top = textGrob("Quebec 2022", 
                                                   vjust = .5, gp = gpar(cex = 1, fontsize = 18)))

qc2022.preds.g

dev.off()

################################################################################
## RANDOM INTERCEPT VS RANDOM SLOPE
################################################################################

anova(qc2022.ri.fit, qc2022.rs.fit.bob) #put the simpler model first

################################################################################
## MARGINAL EFFECTS
################################################################################

qc2022.rs.univ.pred.prop <- ggeffect(qc2022.rs.fit.bob, terms = c("univ"))
qc2022.rs.vote.pred.prop <- ggeffect(qc2022.rs.fit.bob, terms = c("vote_district"))
qc2022.rs.margin.pred.prop <- ggeffect(qc2022.rs.fit.bob, terms = c("margin"))
qc2022.rs.reelected.pred.prop <- ggeffect(qc2022.rs.fit.bob, terms = c("reelected"))
qc2022.rs.int.pred.prop <- ggeffect(qc2022.rs.fit.bob, terms = c("univ", "vote_district"))

qc2022.rs.univ.pred.ref <- ggpredict(qc2022.rs.fit.bob, terms = c("univ"))
qc2022.rs.vote.pred.ref <- ggpredict(qc2022.rs.fit.bob, terms = c("vote_district"))
qc2022.rs.margin.pred.ref <- ggpredict(qc2022.rs.fit.bob, terms = c("margin"))
qc2022.rs.reelected.pred.ref <- ggpredict(qc2022.rs.fit.bob, terms = c("reelected"))
qc2022.rs.int.pred.ref <- ggpredict(qc2022.rs.fit.bob, terms = c("univ", "vote_district"))

ggeffect(qc2022.rs.fit.bob, terms = c("vote_district", "univ")) %>% plot()
ggpredict(qc2022.rs.fit.bob, terms = c("vote_district", "univ")) %>% plot()